setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/bayes_p")
sel1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49")
sel1=paste0(sel1,".snp.bayes_p.txt")
##########这是total_reads>10的部分
i=1
file=read.table(sel1[i],head=T,sep="\t")
file=file[file$exits=="TRUE",]
if(dim(file)[1]>0){
file$unitID=paste(file[,1],file[,2],file[,3],file[,4],sep=":")
file1=file[file$normal_bayes_pvalue>0.1&file$tumor_bayes_pvalue>0.1,]
final=file1[,c("normal_reads1","normal_reads2","normal_var_freq","tumor_reads1","tumor_reads2","tumor_var_freq","unitID")]
final$tumor_var_freq=as.numeric(gsub("%","",final$tumor_var_freq))/100
final$normal_var_freq=as.numeric(gsub("%","",final$normal_var_freq))/100
final$total_reads=as.numeric(final$normal_reads1)+as.numeric(final$normal_reads2)+as.numeric(final$tumor_reads1)+as.numeric(final$tumor_reads2)
sels=which(final$total_reads>=10)
final=final[sels,]
final=final[,-8]
colnames(final)=c(paste(strsplit(gsub(".snp","",sel1[i]),"_")[[1]][1],c("_reads1","_reads2","_var_freq"),sep=""),paste(strsplit(gsub(".snp.bayes_p.txt","",sel1[i]),"_")[[1]][2],c("_reads1","_reads2","_var_freq"),sep=""),"unitID")
}

for (i in 2:length(sel1)){
file=read.table(sel1[i],head=T,sep="\t")
file=file[file$exits=="TRUE",]
if(dim(file)[1]>0){
file$unitID=paste(file[,1],file[,2],file[,3],file[,4],sep=":")
file1=file[file$normal_bayes_pvalue>0.1&file$tumor_bayes_pvalue>0.1,]
final1=file1[,c("normal_reads1","normal_reads2","normal_var_freq","tumor_reads1","tumor_reads2","tumor_var_freq","unitID")]
final1$tumor_var_freq=as.numeric(gsub("%","",final1$tumor_var_freq))/100
final1$normal_var_freq=as.numeric(gsub("%","",final1$normal_var_freq))/100
final1$total_reads=as.numeric(final1$normal_reads1)+as.numeric(final1$normal_reads2)+as.numeric(final1$tumor_reads1)+as.numeric(final1$tumor_reads2)
sels=which(final1$total_reads>=10)
final1=final1[sels,]
final1=final1[,-8]
colnames(final1)=c(paste(strsplit(gsub(".snp","",sel1[i]),"_")[[1]][1],c("_reads1","_reads2","_var_freq"),sep=""),paste(strsplit(gsub(".snp.bayes_p.txt","",sel1[i]),"_")[[1]][2],c("_reads1","_reads2","_var_freq"),sep=""),"unitID")

final=merge(final,final1,by="unitID",all=T)
}
}
length(unique(final$unitID))
[1] 12284

write.table(final,"nobias_AShM_total_reads10.txt",quote=F,row.names = F,sep="\t")

##########这是不对total_reads筛选的部分
i=1
file=read.table(sel1[i],head=T,sep="\t")
file=file[file$exits=="TRUE",]
if(dim(file)[1]>0){
file$unitID=paste(file[,1],file[,2],file[,3],file[,4],sep=":")
file1=file[file$normal_bayes_pvalue>0.1&file$tumor_bayes_pvalue>0.1,]
final=file1[,c("normal_reads1","normal_reads2","normal_var_freq","tumor_reads1","tumor_reads2","tumor_var_freq","unitID")]
final$tumor_var_freq=as.numeric(gsub("%","",final$tumor_var_freq))/100
final$normal_var_freq=as.numeric(gsub("%","",final$normal_var_freq))/100
final$total_reads=as.numeric(final$normal_reads1)+as.numeric(final$normal_reads2)+as.numeric(final$tumor_reads1)+as.numeric(final$tumor_reads2)
#sels=which(final$total_reads>=10)
#final=final[sels,]
final=final[,-8]
colnames(final)=c(paste(strsplit(gsub(".snp","",sel1[i]),"_")[[1]][1],c("_reads1","_reads2","_var_freq"),sep=""),paste(strsplit(gsub(".snp.bayes_p.txt","",sel1[i]),"_")[[1]][2],c("_reads1","_reads2","_var_freq"),sep=""),"unitID")
}

for (i in 2:length(sel1)){
file=read.table(sel1[i],head=T,sep="\t")
file=file[file$exits=="TRUE",]
if(dim(file)[1]>0){
file$unitID=paste(file[,1],file[,2],file[,3],file[,4],sep=":")
file1=file[file$normal_bayes_pvalue>0.1&file$tumor_bayes_pvalue>0.1,]
final1=file1[,c("normal_reads1","normal_reads2","normal_var_freq","tumor_reads1","tumor_reads2","tumor_var_freq","unitID")]
final1$tumor_var_freq=as.numeric(gsub("%","",final1$tumor_var_freq))/100
final1$normal_var_freq=as.numeric(gsub("%","",final1$normal_var_freq))/100
final1$total_reads=as.numeric(final1$normal_reads1)+as.numeric(final1$normal_reads2)+as.numeric(final1$tumor_reads1)+as.numeric(final1$tumor_reads2)
#sels=which(final1$total_reads>=10)
#final1=final1[sels,]
final1=final1[,-8]
colnames(final1)=c(paste(strsplit(gsub(".snp","",sel1[i]),"_")[[1]][1],c("_reads1","_reads2","_var_freq"),sep=""),paste(strsplit(gsub(".snp.bayes_p.txt","",sel1[i]),"_")[[1]][2],c("_reads1","_reads2","_var_freq"),sep=""),"unitID")

final=merge(final,final1,by="unitID",all=T)
}
}
length(unique(final$unitID))
[1] 21421

write.table(final,"nobias_AShM.txt",quote=F,row.names = F,sep="\t")

